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Abstract 

This paper develops a novel framework for phase retrieval, a problem which arises in X-ray crystal- 
lography, diffraction imaging, astronomical imaging and many other applications. Our approach, called 
PhaseLift, combines multiple structured illuminations together with ideas from convex programming to 
recover the phase from intensity measurements, typically from the modulus of the diffracted wave. We 
demonstrate empirically that any complex-valued object can be recovered from the knowledge of the 
magnitude of just a few diffracted patterns by solving a simple convex optimization problem inspired by 
the recent literature on matrix completion. More importantly, we also demonstrate that our noise-aware 
algorithms are stable in the sense that the reconstruction degrades gracefully as the signal-to-noise ratio 
decreases. Finally, we introduce some theory showing that one can design very simple structured illumi- 
nation patterns such that three diffracted figures uniquely determine the phase of the object we wish to 
recover. 

Keywords. Diffraction, Fourier transform, convex optimization, trace-norm minimization. 



1 Introduction 

1.1 The phase retrieval problem 

This paper considers the fundamental problem of recovering a general signal, an image for example, from the 
magnitude of its Fourier transform. This problem, also known as phase retrieval, arises in many applications 
and has challenged engineers, physicists, and mathematicians for decades. Its origin comes from the fact that 
detectors can often times only record the squared modulus of the Fresnel or Fraunhofer diffraction pattern 
of the radiation that is scattered from an object. In such settings, one cannot measure the phase of the optical 
wave reaching the detector and, therefore, much information about the scattered object or the optical field is 
lost since, as is well known, the phase encodes a lot of the structural content of the image we wish to form. 



Historically, the first application of phase retrieval is X-ray crystallography 1 32 52 1, and today this may 
still very well be the most important application. Over the last century or so, this field has developed a wide 
array of techniques to recover Bragg peaks from missing-phase data. Of course, the phase retrieval problem 
permeates many other areas of imaging science, and other applications include diffraction imaging |jT2l, 



optics |65 1, astronomical imaging 1 20 1, microscopy |49 1, to name just a few. In particular, X-ray tomography 



has become an invaluable tool in biomedical imaging to generate quantitative 3D density maps of extended 
specimens on the nanoscale [21]. Other subjects where phase retrieval plays an important role are quantum 
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mechanics 1 19 59 1 and even differential geometry |[8|. We note that phase retrieval has seen a resurgence in 



activity in recent years, fueled on the one hand by the desire to image individual molecules and other nano- 
particles, and on the other hand by new imaging capabilities: one such recent modality is the availability of 
new X-ray synchrotron sources that provide extraordinary X-ray fluxes, see for example |[9j|2Tj49[53j62|. 



References and various instances of the phase retrieval problem as well as some theoretical and numerical 
solutions can be found in |35[|40[[44j . 



There are many ways in which one can pose the phase-retrieval problem, for instance depending upon 
whether one assumes a continuous or discrete-space model for the signal. In this paper, we consider finite 
length signals (one-dimensional or multi-dimensional) for simplicity, and because numerical algorithms 
ultimately operate with digital data. To fix ideas, suppose we have a ID signal x = (x[0], a:[l], . . . ,x[n — 
1] ) G C" and write its Fourier transform as 



X UJ 



^ x[t]e-^2-^*/", uen. (1.1) 

0<t<n 



Here, $7 is a grid of sampled frequencies and an important special case is = {0, 1, . . . , n — 1} so that 
the mapping is the classical unitary discrete Fourier transform (DFT|^ The phase retrieval problem consists 
in finding x from the magnitude coefficients a; G 0,. When Q, is the usual frequency grid as above 

and without further information about the unknown signal x, this problem is ill-posed since there are many 
different signals whose Fourier transforms have the same magnitude. Clearly, if x is a solution to the phase 
retrieval problem, then (1) cx for any scalar c G C obeying |c| = 1 is also solution, (2) the "mirror function" 
or time-reversed signal x[—t mod n] where t = 0, 1, . . . , n — 1 is also solution, and (3) the shifted signal 
x[t — a mod n] is also a solution. From a physical viewpoint these "trivial associates" of x are acceptable 
ambiguities. But in general infinitely many solutions can be obtained from {|x[a;]| : w G il} beyond these 



trivial associates |61 1. 



1.2 Main approaches to phase retrieval 

Holographic techniques are among the more popular methods that have been proposed to measure the phase 
of the optical wave. While holographic techniques have been successfully applied in certain areas of optical 
imaging, they are generally difficult to implement in practice p2| . Hence, the development of algorithms for 
signal recovery from magnitude measurements is still a very active field of research. Existing methods for 
phase retrieval rely on all kinds of a priori information about the signal, such as positivity, atomicity, support 
constraints, real-valuedness, and so on [18| [26l[27l|47| . Direct methods [33 1 are limited in their applicability 
to small-scale problems due to their large computational complexity. 

Oversampling in the Fourier domain has been proposed as a means to mitigate the non-uniqueness 
of the phase retrieval problem. While oversampling offers no benefit for most one-dimensional signals, the 
situation is more favorable for multidimensional signals, where it has been shown that twofold oversampling 
in each dimension almost always yields uniqueness for finitely supported, real-valued and non-negative 



signals |11 34 61|. In other words, a digital image of the form x = {x[ti,t2]} with < ti < ni and 



< t2 < n2, whose Fourier transform is given by 



x[uji,uj2] = — ^ y 2;[ti,t2]e-'^^"('^i*i/"^+'^2*^/"2\ (1.2) 



'For later reference, we denote the Fourier transform operator by F and the inverse Fourier transform by F ^ . 
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is usually uniquely determined from the values of |x[a;i,a;2]| on the oversampled grid uj = {uji,uj2) G = 
r^i X 0.2 in which Q.i = {0, 1/2, 1, 3/2, . . . , + 1/2}. This holds provided x has proper spatial support, is 
real valued and non-negative. 

As pointed out in [44], these uniqueness results do not say anything about how a signal can be recovered 
from its intensity measurements, or about the robustness and stability of commonly used reconstruction 
algorithms — a fact we shall make very clear in the sequel. In fact, theoretical uniqueness conditions do not 
readily translate into numerical methods and as a result, the algorithmical and practical aspects of the phase 
retrieval problem (from noisy data) still pose significant challenges. 

By and large, the most popular methods for phase retrieval from oversampled data are alternating projec- 



tion algorithms pioneered by Gerchberg and Saxton [29] and Fienup p6] 27 1. These methods often require 
careful exploitation of signal constraints and delicate parameter selection to increase the likelihood of con- 
vergence to a correct solution [ 18,46,47 ,56]. We describe the simplest realization of a widely used approach 




based on alternating projections |50 1, which assumes support constraints in the spatial domain and oversam- 
pled measurements in the frequency domain. With T being a known subset containing the support of the 
signal X (supp(x) C T) and Fourier magnitude measurements {y[oj\]tjj^Q^ with y[uj\ = |x[(x)]|, the method 
works as follows: 

1. Initialization: Choose an initial guess xq and set zqIu:] = y[uj\ for u; G 17. 

2. Loop: For A; = 1, 2, . . . inductively define 

(1) Xk[t] = 



(2) ZfcH=yM^^ fovuen 
until convergence. 

While this algorithm is simple to implement and amenable to additional constraints such as the positivity 
of X, its convergence remains problematic. Projection algorithms onto convex sets are well understood 
|[3j[Toj|3T]|66[. However, the set {z : \z[uj]\ = \x[uj]\} is not convex and, therefore, the algorithm is not 
known to converge in general or even to give a reasonable solution |[3[|4l]|44|. Good results have been 
reported in certain settings but they appear to be nevertheless somewhat problematic in light of our numerical 
experiments from Section]?] Moreover, as discussed in [38] , one of the most stringent limitations of these 
methods is the need for isolated objects (the support constraint). Finally, | ]45| points out that oversampling 
is not always practically feasible as certain experimental geometries allow only for sub-Nyquist sampling; 
an example is the Bragg sampling from periodic crystalline structures. 

In a different direction, a frame-theoretic approach to phase retrieval has been proposed in |JT]]2|, where 
the authors derive various necessary and sufficient conditions for the uniqueness of the solution, as well as 
various numerical algorithms. While theoretically appealing, the practical applicability of these results is 
limited by the fact that very specific types of measurements are required, which cannot be realized in most 
applications of interest. 

To summarize our discussion, we have seen many methods which all represent some important attempts 
to find efficient algorithms, and work well in certain situations. However, these techniques do not always 
provide a consistent and robust result. 
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1.3 PhaseLift - a novel methodology 



This paper develops a novel methodology for phase retrieval based on a rigorous and flexible numerical 
framework. Whereas most of the existing methods seek to overcome nonuniqueness by imposing additional 
constraints on the signal, we pursue a different direction by assuming no constraints at all on the signal. 
There are two main components to our approach. 

• Multiple structured illuminations. We suggest collecting several diffraction patterns providing 'differ- 
ent views' of the sample or specimen. This can be accomplished in a number of ways: for instance, 
by modulating the light beam falling onto the sample or by placing a mask right after the sample, see 
Section [2] for details. Taking multiple diffraction patterns usually yields uniqueness as discussed in 
Section [3l 

The concept of using multiple measurements as an attempt to resolve the phase ambiguity for diffrac- 
tion imaging is of course not new, and was suggested in [54|. Since then, a variety of methods have 
been proposed to carry out these multiple measurements; depending on the particular application, 
these may include the use of various gratings and/or of masks, the rotation of the axial position of the 
sample, and the use of defocusing implemented in a spatial light modulator, see [22] for details and 
references. Other approaches include ptychography, an exciting field of research, where one records 
several diffraction patterns from overlapping areas of the sample, see ||60 63 1 and references therein. 



• Formulation of phase recovery as a matrix completion problem. We suggest (1) lifting up the prob- 
lem of recovering a vector from quadratic constraints into that of a recovering of a rank-one matrix 
from affine constraints, and (2) relaxing the combinatorial problem into a convenient convex program. 
Since the lifting step is fundamental to our approach, we will refer to the proposed numerical frame- 
work as PhaseLift. The price we pay for trading the nonconvex quadratic constraints into convex 
constraints is that we must deal with a highly underdetermined problem. However, recent advances in 
the areas of compressive sensing and matrix completion have shown that such convex approximations 
are often exact. 

Although our algorithmic framework appears to be novel for phase retrieval, the idea of solving prob- 
lems involving nonconvex quadratic constraints by semidefinite relaxations has a long history in opti- 
mization, see ||7| and references therein, and Section [T!4l below for more discussion. 

The goal of this paper is to demonstrate that taken together, multiple coded illuminations and convex pro- 
gramming (trace-norm minimization) provide a powerful new approach to phase retrieval. Further, a signif- 
icant aspect of our methodology is that our systematic optimization framework offers a principled way of 
dealing with noise, and makes it easy to handle various statistical noise models. This is important because 
in practice, measurements are always noisy. In fact, our framework can be understood as an elaborate regu- 
larized maximum likelihood method. Lastly, our framework can also include a priori knowledge about the 
signal that can be formulated or relaxed as convex constraints. 



1.4 Precedents 

At the abstract level, the phase retrieval problem is that of finding x E C" obeying quadratic equations of the 
form \ {ak, x)p = 6^- Casting such quadratic constraints as affine constraints about the matrix variable X = 
XX* has been widely used in the optimization literature for finding good bounds on a number of quadratically 
constrained quadratic problems (QCQP). Indeed, solving the general case of a QCQP is known to be an NP- 
hard problem since it includes the family of boolean linear programs [7]. The approach usually consists 
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in finding a relaxation of the QCQP using semidefinite programming (SDP), for instance via Lagrangian 
duality. An important example of this strategy is Max Cut, an NP-hard problem in graph theory which can 
be formulated as a QCQP In a celebrated paper, Goemans and Williamson introduced a relaxation |30| for 
this problem, which lifts or linearizes a nonlinear, nonconvex problem to the space of symmetric matrices. 
Although there are evident connections to our work, our relaxation is quite different from these now standard 
techniques. 

The idea of linearizing the phase retrieval problem by reformulating it as a problem of recovering a 
matrix from linear measurements can be found in ||T|. While this reference also contains some intriguing 
numerical recovery algorithms, their practical relevance for most applications is limited by the fact that the 
proposed measurement matrices either require a very specific algebraic structure which does not seem to be 
compatible with the physical properties of diffraction, or the number of measurements is proportional to the 
square of the signal dimension, which is not feasible in most applications. 

In terms of framework, the closest approach is the paper 1 17 1, in which the authors use a matrix com- 
pletion approach for array imaging from intensity measurements. Although this paper executes a similar 
relaxation as ours, there are some differences. We present a "noise-aware" framework, which makes it pos- 
sible to account for a variety of noise models in a systematic way. Moreover, our emphasis is on a novel 
combination of structured illuminations and convex programming, which seems to bear great potential. 



2 Methodology 

2.1 Structured illumination 

Suppose X = {x[t]} is the object of interest (t may be a one- or multi-dimensional index). In this paper, 
we shall discuss illumination schemes collecting the diffraction pattern of the modulated object 
where the waveforms or patterns w[t] may be selected by the user. There are many ways in which this can 
be implemented in practice, and we discuss just a few of those. 



Masking. One possibility is to modify the phase front after the sample by inserting a mask or a phase 
plate, see |42| for example. A schematic layout is shown in Figure[T] In |38 1, the sample is scanned 



by shifting the phase plate as in ptychography (discussed below); the difference is that one scans the 
known phase plate rather than the object being imaged. 

Optical grating. Another standard approach would be to change the profile or modulate the illumi- 
nating beam, which can easily be accomplished by the use of optical gratings [43 1. A simplified 
representation would look similar to the scheme depicted in Figure [TJ with a grating instead of the 
mask (the grating could be placed before or after the sample). 

Ptychography. Here, one measures multiple diffraction patterns by scanning a finite illumination on 



an extended specimen 1 60 63 1 . In this setup, it is common to maintain a substantial overlap between 



adjacent illumination positions. 

Oblique illuminations. One can use illuminating beams hitting the sample at user specified angle | [23| , 
see Figure [2] for a schematic illustration of this approach. One can also imagine having multiple 
simultaneous oblique illuminations. 



As is clear, there is no shortage of options and one might prefer solutions which require generating as 
few diffraction patterns as possible for stable recovery. 
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Figure 1: A typical setup for structured illuminations in diffraction imaging using a phase mask. 




Figure 2: A typical setup for structured illuminations in diffraction imaging using oblique illumina- 
tions. The left image shows direct (on-axis) illumination and the right image corresponds to obhque 
(off-axis) illumination. 
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2.2 Lifting 

Suppose we have xq G C" or C"^^"^ (or some higher-dimensional version) about which we have quadratic 
measurements of the form 

A(xo) = {Kafc,Xo)|^ : /c = l,2,...,m}. (2.1) 

In the setting where we would collect the diffraction pattern of zi;[t]xo[t] as discussed earlier, then the 
waveform a]^[t] can be written as 



here, is a frequency value so that au [t] is a patterned complex sinusoid. One can assume for convenience 
that the normalizing constant is such that is unit normed, i.e. \\ak\W = l^feMP = 1- Phase retrieval is 
then the feasibility problem 

^"'^ ^ (2 2) 

obeying A(x) = A(xo) := h. 

As is well known, quadratic measurements can be lifted up and interpreted as linear measurements about 
the rank-one matrix X = xx* . Indeed, 

|(afc,x)p = Ti{x*akalx) = Ti{akalxx*) :=Tr(^fe^), 

where Ak is the rank-one matrix a^a^. In what follows, we will let A be the linear operator mapping positive 
semidefinite matrices into {Tr{AkX) : A: = 1, . . . , m}. Hence, the phase retrieval problem is equivalent to 



find X 

subject to A{X) = i 
X ^ 
rank(X) 



minimize rank(X) 
4» subject to A{X) = h (2.3) 
X ^ 0. 

Upon solving the left-hand side of ( |2.3| ), we would factorize the rank-one solution X as xx* , hence finding 
solutions to the phase-retrieval problem. Note that the equivalence between the left- and right-hand side of 



( |2.3| ) is straightforward since by definition, there exists a rank-one solution. Therefore, our problem is a rank 
minimization problem over an affine slice of the positive semidefinite cone. As such, it falls in the realm of 
low-rank matrix completion or matrix recovery, a class of optimization problems that has gained tremendous 



attention in recent years, see e.g. 1 13 14 58 1. Just as in matrix completion, the linear system A{X) = b, 
with unknown in the positive semidefinite cone, is highly underdetermined. For instance suppose our signal 
xq has n complex unknowns. Then we may imagine collecting six diffraction patterns with n measurements 
for each (no oversampling). Thus m = 6n whereas the dimension of the space of n x n Hermitian matrices 
over the reals is n^, which is obviously much larger. 

We are of course interested in low-rank solutions and this makes the search feasible. This also raises an 
important question: what is the minimal number of diffraction patterns needed to recover x, whatever x may 
be? Since each pattern yields n real-valued coefficients and x has n complex-valued unknowns, the answer 
is at least two. Further, in the context of quantum state tomography. Theorem II in 1^28] shows one needs 
at least 3n — 2 intensity measurements to guarantee uniqueness, hence suggesting an absolute minimum of 
three diffraction patterns. Are three patterns sufficient? For some answers to this question, see Section|3] 
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2.3 Recovery via convex programming 



The rank minimization problem (2.3) is NP hard. We propose using the trace norm as a convex surrogate 



for the rank functional, giving the familiar SDP (and a crucial component of PhaseLift), 

minimize trace(X) 

subject to A{X) = b (2.4) 

X to. 

This problem is convex and there exists a wide array of numerical solvers including the popular Nesterov's 



accelerated first order method p5|. As far as the relationship between (2.3 1 and (|2.4[) is concerned, the 



matrix A in most diffraction imaging applications is not known to obey any of the conditions derived in the 



literature 1 13 14 58| that would guarantee a formal equivalence between the two programs. Nevertheless, 



the formulation ( 2.4| ) enjoys great empirical performance as demonstrated in Section[4] 



We mentioned earlier that measurements are typically noisy and that our formulation allows for a prin- 
cipled approach to deal with this issue for a variety of noise models. Suppose the measurement vector 
{bk} is sampled from a probability distribution p{-; n), where /i = A(xo) is the vector of noiseless values, 

= Kq^/c) 2;o)P- Then a classical fitting approach simply consists of maximizing the Ukelihood, 

maximize p{b\lj) 
subject to /X = A(x) 

with optimization variables /i and x. (A more concise description is to find x such that p{b] A(x)) is maxi- 
mum.) This is, of course, not tractable and our convex formulation suggests solving instead 

minimize — logp(6; ^) + A Tr(X) 

subject to 11 = A{X) (2.6) 
X ^ 

with optimization variables ^ and X (in other words, find X ^ such that — logp(6; A{X)) + A Tr(X) is 
minimum). Above, A is a positive scalar and, hence, our approach is a penalized or regularized maximum 
likelihood method, which trades off between goodness and complexity of the fit. When the likelihood is 
log-concave, problem ( |2.6[ ) is convex and solvable. We give two examples for concreteness: 

• Poisson data. Suppose that {bk} is a sequence of independent samples from the Poisson distributions 
Poi(^fc). The Poisson log-likelihood for independent samples has the form bk log /x^ — ^k (up to 
an additive constant factor) and thus, our problem becomes 

minimize Y.k \^^k - bk log Pk] + A Tr (X) 
subject to /i = A{X) 

X yo. 

• Gaussian data. Suppose that {6^} is a sequence of independent samples from the Gaussian distribu- 
tion with mean fik and variance (t| (or is well approximated by Gaussian variables). Then our problem 
becomes 

minimize ^ (bk - fJ-k)"^ + A Tr{X) 

k 

subject to /i = A{X) 
X ^ 0. 
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If S is a diagonal matrix with diagonal elements o"^, this can be written as 

minimize l[b - A{X)]*^''^[b - A{X)] + XTt{X) 
subject to X 0. 

Both formulations are of course convex and in both cases, one recovers the noiseless trace-minimization 



problem ( |X4| as A 0+. 

In addition, it is straightforward to include further constraints frequently discussed in the phase retrieval 
literature such as real-valuedness, positivity, atomicity and so on. Suppose the support of x is known to be 
included in a set T known a priori. Then we would add the linear constraint 

Xij = 0, {i,j)^TxT. 

(Algorithmically, one would simply work with a reduced-size matrix.) Suppose we would like to enforce 
real-valuedness, then we simply assume that X is real valued and positive semidefinite. Finally positivity 
can be expressed as linear inequalities 

Xij > 0. 

Of course, many other types of constraints can be incorporated in this framework, which provides apprecia- 
ble flexibility. 

2.4 PhaseLift with reweighting 

The trace norm promotes low-rank solutions and this is why it is often used as a convex proxy for the 
rank. However, it is possible to further promote low-rank solutions by solving a sequence of weighted trace- 
norm problems, a technique which has been shown to provide even more accurate solutions |15]25|. The 



reweighting scheme works like this: choose e > 0; start with Wq = I and for A; = 0, 1, ... , inductively 
define Xk as the optimal solution to 

minimize trace (VF^X) 

subject to A{X) = b (2.7) 
X ^ 



and update the 'weight matrix' as 



Wk+i = {Xk + eI)-\ 



The algorithm terminates on convergence or when the iteration count k attains a specified maximum number 
of iterations k^ax- One can see that the first step of this procedure is precisely ( |2.4[ ); after this initial step, 
the algorithm proceeds in solving a sequence of trace-norm problems in which the matrix weights Wk are 
roughly the inverse of the current guess. 



As explained in the literature |24 25 1, this reweighting scheme can be viewed as attempting to solve 



minimize f{X) = log(det(X + el)) 

subject to A{X) = b (2.8) 
X ^ 



by minimizing the tangent approximation to / at each iterate; that is to say, at step k, ( |2.7| ) is equivalent to 
minimizing f{Xk-i) + {V f{Xf^_i), X — X^-i) over the feasible set. This can also be applied to noise- 



aware variants where one would simply replace the objective functional in ( |2.6[ ) with 

-logp{b;fi) + XTT{WkX), 
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at each step, and update Wk in exactly the same way as before. 

The log-det functional is closer to the rank functional than the trace norm. In fact, minimizing this 
functional solves the phase retrieval problem as incorporated in the following theorem. 

Theorem 2.1 Suppose A is one to one and that the identity matrix I is in the span of the sensing matrices 



Aj.. Then the unique solution of the phase retrieval problem ( |2.2[ ) is also the unique minimizer to (2.8 1 — up 
to global phase^ This holds for all values of £ > 0. The same conclusion holds without the inclusion 
assumption provided one modifies the reweighting scheme and substitutes the objective function f{X) in 
g with f {RXR) where R = (^^ AkY'^. 



Since the reweighting algorithm is a good heuristic for solving ( |2.8| ), we potentially have an interesting and 
tractable method for phase retrieval. It is not a perfect heuristic, however, as we cannot expect this procedure 
to always find the global minimum since the objective functional is concave. 

The assumption that the identity matrix is in the span of the A^'s holds whenever the modulus of the 
Fourier transform of the sample is measured. Indeed, if is observed, then letting {/^} be the rows of 
F, we have Ekfkf*k=I- 

Proof Our assumption implies that for any feasible X, Tt{X) = J2k Tr(^fcX) = J2k ^k^k is fixed. 
Assume without loss of generality that feasible points obey Tt{X) = 1 (if Ti{X) = 0, then the unique 
solution is X = 0). If xq is the unique solution to phase retrieval (up to global phase), then Xq = xqXq 
is the only rank-one feasible point. We thus need to show that any feasible X with rank(X) > 1, obeys 
f{X) > /(Xq), a fact which follows from the strong concavity of / (of the logarithm). Let X = J2j ^j'^j'^) 
be any eigenvalue decomposition of a feasible point. Then 

f{X) = log(det(e/ + X)) = Y, \og{e + A,), 

3 

and it follows from the strict concavity of the log that 

^log(e + \j) > Xj log(e + 1) + (1 - Xj) logs = log(e + 1) + (n - 1) lege. 

j 3 

The first strict inequality holds unless X is rank one, in which case, we have equality. The equality follows 
from Y^- Xj = Tr(X) = 1. Since the right-hand side is nothing else than /{Xq), the theorem is established. 

For the second part, set Bk = R^^AkR^^ and consider a new data problem with constraints {X : 
TT{BkX) = bk,X >z 0}. Now X is feasible for our problem if and only if RXR is feasible for this new 
problem. (This is because the mapping X i— )• RXR preserves the positive semidefinite cone.) Now suppose 
Xq is the solution to phase retrieval and set Xq = xqXq as before. Since / is in the span of the sensing 
matrices Bk, we have just learned that for all for all RXR ^ RXqR and X feasible for our problem, 

f{RXR) > f{RXoR). 

This concludes the proof. ■ 

^Quadratic measurements can of course never distinguish between x and cx in which c £ C has unit norm. When the solution 
is unique up to a multiplication with such a scalar, we say that unicity holds up to global phase. From now on, whenever we talk 
about unicity, it is implied up to global phase. 
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3 Theory 



Our PhaseLift framework poses two main theoretical questions: 

1 . When do multiple diffracted images imply unicity of the solution? 

2. When does our convex heuristic succeed in recovering the unique solution to the phase-retrieval prob- 
lem? 

Developing comprehensive answers to these questions constitutes a whole research program, which clearly 
is beyond the scope of this work. In this paper, we shall limit ourselves to introducing some theoretical 
results showing simple ways of designing diffraction patterns, which give unicity. Our focus is on getting 
uniqueness from a very limited number of diffraction patterns. For example, we shall demonstrate that in 
some cases three diffraction images are sufficient for perfect recovery. Thus, we give below partial answers 
to the first question and will address the second in a later publication. 

A frequently discussed approach to retrieve phase information uses a technique from holography. Roughly 
speaking, the idea is to let the signal of interest x interfere with a known reference beam y. One typically 
measures \x + and \x — iy\'^ and precise knowledge of y allows, in principle, to recover the amplitude 
and phase of x. Holographic techniques are hard to implement \22\ in practice. Instead, we propose using a 
modulated version of the signal itself as a reference beam which in some cases may be easier to implement. 

To discuss this idea, we need to introduce some notation. For a complex signal z G C", we let |zp 
be the nonnegative real-valued n-dimensional vector containing the squared magnitudes of z. Suppose first 
that X is a one-dimensional signal (x[0], . . . , x[n — 1]) and that F„ is the n x n unitary DFT. In this 
section, we consider taking 3n real- valued measurements of the form 

A(x) = \Fnix + D'x)\^, \Fn{x - iD'x)\^}, (3.1) 

where D is the modulation 

Z? = diag({e*2-*/"}o<t<„-i). 

These measurements can be obtained by illuminating the sample with the three light fields 1, e*^'^^*/" and 
gj27r(st/n-i/4) show bclow that thcsc 3n measurements are generally sufficient for perfect recovery. 

Theorem 3.1 Suppose the DFT of x £ C" does not vanish. Then x can be recovered up to global phase 
from the 3n real numbers A(x) ( |3.1| ) if and only if s is prime with n. In particular, assuming primality, if the 
trace-minimization program ( |2.4| ) or the iteratively reweighted algorithm return a rank-1 solution, then this 
solution is exact. 

Conversely, if the DFT vanishes at two frequency points k and k' obeying k — k' ^ s mod n, then 
recovery is not possible, from the 3n real numbers p.l[ ). 

The proof of this theorem is constructive and we give a simple algorithm that achieves perfect reconstruction. 
Further, one can use masks to scramble the Fourier transform as to make sure it does not vanish. Suppose 
for instance that we collect 

A{Wx), W = dmg{{z[t]}o<t<n-i). 

where the z[t]'s are iid A/'(0, 1). Then since the Fourier transform of does not vanish with probability 

one, we have the following corollary. 

Corollary 3.2 Assume s is prime with n. Then with probability one, x can be recovered up to global phase 
from the 3n real numbers K{Wx) where W is the diagonal matrix with Gaussian entries above. 
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Of course, one could derive similar results by scrambling the Fourier transform with the aid of other types 
of masks, e.g. binary masks. We do not pursue such calculations. 

We now turn our attention to the situation in higher dimensions and will consider the 2D case (higher 
dimensions are treated in the same way). Here, we have a discrete signal ^2] £ c^ix^a about which 
we take the 3nin2 measurements 

{\FniXn2x\'^ ■, \J^niXn2{x + 'D''x)\^, \FniXn2{x - W x)\^] , S = (si,S2); (3.2) 

J^nixn2 is the 2D unitary Fourier transform defined by ( |1.2[ ) in which the frequencies belong to the 2D grid 
{0, 1, . . . , ni — 1} X {0, 1, . . . , n2 — 1}, and is the modulation 

With these definitions, we have the following result: 

Theorem 3.3 Suppose the DFTofx G (7"i x"2 does not vanish. Then x can be recovered up to global phase 
from the 3nin2 real numbers (|3.2|) if and only if si is prime with ni, S2 is prime with n2 and rii is prime 



with n2. Under these assumptions, if the trace-minimization program (2.4 1 or the iteratively reweighted 
algorithm return a rank-1 solution, then this solution is exact. 

Again, one can apply a random mask to turn this statement into a probabilistic statement holding either 
with probability one or with very large probability depending upon the mask that is used. 

One can always choose si and S2 such that they be prime with ni and n2 respectively. The last condition 
may be less friendly but one can decide to pad one dimension with zeros to guarantee primality. This 
is equivalent to a slight oversampling of the DFT along one direction. An alternative is to take 5nin2 
measurements in which we modulate the signal horizontally and then vertically; that is to say, we modulate 
with s = (si, 0) and then with s = (0, S2). These 5nin2 measurements guarantee recovery if si is prime 



with ni and S2 is prime with n2 for all sizes ni and n2, see Section 3.3 for details. 



3.1 Proof of Theorem O 

Let X = (£[0], . . . , x[n — 1]) be the DFT of x. Then knowledge of A(2;) is equivalent to knowledge of 

\x[k] + x[k - s]p, and \x[k\ - ix[k - s]p 

for all A; G {0, 1, . . . , n — 1} (above, k — s is understood mod n). Write x[k] = |x[A:]|e*'^['^l so that (()[k] is 
the missing phase, and observe that 

\x[k] + x[k - s]|2 = \x[k]\^ + \x[k - s]p + 2\x[k]\\x[k - s]|Re(e*(<^['=-^]-<^W) 
\x[k] - ix[k - s]|2 = |x[A;]|2 + \x[k - s]|2 + 2\x[k]\\x[k - s]|Im(e^('^t^-^l-<^W). 

Hence, if x[k] / for all A; € {0, 1, . . . , n — 1}, our data gives us knowledge of all phase shifts of the form 

(t)[k - s] - 4>[k], k = 0,l,...,n-l. 

We can, therefore, initialize (/}[()] to be zero and then get the values of 0[— s], 2s] and so on. 

This process can be represented as a cycle in the group Z/nZ as the sequence (0, — s, —2s, . . .). We 
would like this cycle to contain n unique elements, which is true if and only if the cyclic subgroup (0, s, 2s, . . .) 
has order n. This is equivalent to requiring gcd(s, n) = 1. If this subgroup has a smaller order, then recovery 
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is impossible since we finish the cycle before we have all the phases; the phases that we are able to recover 
do not enable us to determine any more phases without making further assumptions. 

For the second part of the theorem, assume without loss of generality, that s = — 1 and that {k,k') = 
(0, ko) (1 < ko < n — 1). For simplicity suppose these are the only zeros of the DFT. This creates two 
disjoint sets of frequency indices: those for which < k < ko and those for which ko < k < n — I. We 
are given no information about the phase difference between elements of these two subgroups, and hence 
recovery is not possible. This argument extends to situations where the DFT vanishes more often, in which 
case, we have even more indeterminacy. 

3.2 Proof of Theorem |33] 

Let X = {x[ki, k2]}, where {ki, k2) € {0, 1, ... ni — 1} x {0, 1, . . . , n2 — 1} be the DFT of x. Then we 
have knowledge of 

\x[ki, k2]\'^, \x[ki,k2] + x[ki — Ss, k2 — S2]|^, and \x[ki, k2] — ix[ki — Si,k2 — S2]P 
for all {kik2). With the same notations as before, this gives knowledge of all phase shifts of the form 

4>[ki - si, k2 - S2] - 4>[ki, k2], < A;i < ni, < ^2 < ^2 - 1. 

Hence, we can initialize </>[0,0] to be zero and then get the values of 0[— si,— S2], (/>[— 2si, — 2s2] and 
so on. The argument is as before: we would like the cyclic subgroup ((0, 0), (si, S2), (2si, 2s2), . . .) in 
Z/niZ X Z/n2Z to have order nin2. Now the order of an element (si, S2) G 'L/ni'L x Z/n2Z is equal to 

lcm( I si I , I S2 1 ) = lcm(ni /gcd(ni , si ) , n2/gcd(n2 ,52)), 

where is the order of si in Z/niZ and likewise for |s2|- Noting that lcm(a, b) < ah and that equality is 
achieved if and only if gcd(a, b) = 1, we must simultaneously have 

gcd(si,ni) = 1, gcd(s2,n2) = l and gcd(ni, 712)) = 1 

to have uniqueness. 

3.3 Extensions 

It is clear from our analysis that if we were to collect \Fniy.n2x\^ together with 

{\J'n^xn2{x + V"'X)\^, |-F„ix„2(x - k = l,...,K, 

SO that one collects {2K + l)nin2 measurements, then 2D recovery is possible if and only if {si, . . . , sk} 
generates Z/niZ x Z/n2Z (and the Fourier transform has no nonzero components). This can be understood 
by analyzing the generators of the group Z/niZ x Z/n2Z. 

A simple instance consists in choosing one modulation pattern to be (si, 0) and another to be (0, 82)- 
If si is prime is ni and S2 with n2, these two modulations generate the whole group regardless of the 
relationship between ni and n2. An algorithmic way to see this is as follows. Initialize </)(0,0). Then 
by using horizontal modulations, one recovers all phases of the form (/)(fci,0). Further, by using vertical 
modulations (starting with 4){ki,Q)), one can recover all phases of the form (l){ki, ^2) by moving upward. 
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4 Numerical Experiments 



This section introduces numerical simulations to illustrate and study the effectiveness of PhaseLift. 



4.1 Numerical solvers 

All numerical algorithms were implemented in Matlab using TFOCS ||6| as well as modifications of TFOCS 
template files. TFOCS is a library of Matlab-files designed to facilitate the construction of first-order meth- 
ods for a variety of convex optimization problems, which include those we consider. 
In a nutshell, suppose we wish to solve the problem 

minimize g{X) := -i{b;A{X)) + ATr(X) 
subject to X ^ 

in which i{b;A{X)) is a smooth and concave (in X) log-likelihood. Then a projected gradient method 
would start with an initial guess Xq, and inductively define 

Xk = V{Xk-i-tkVg{Xk-i)), 

where {tk} is a sequence of stepsize rules and V is the projection onto the positive semidefinite cone. 
(Various stepsize rules are typically considered including fixed stepsizes, backtracking line search, exact 
line search and so on.) 

TFOCS implements a variety of accelerated first-ordered methods pioneered by Nesterov, see p5j and 
references therein. One variant [4] works as follows. Choose Xq, set Yq = Xq and 6*0 = 1, and inductively 
define 

Xk = V{Yk-i - tkVg{Yk-i)) 

A = 0k{Ok'i - 1) 

Yk = + PkiXk — Xk^i) 

where {t^} is a sequence of stepsize rules as before. The sequence {9^} is usually referred to as a sequence 
of accelerated parameters, and {Ifc} is an auxiliary sequence at which the gradient is to be evaluated. The 
advantage of this approach is that the computational work per iteration is as in the projected gradient method 
but the number of iterations needed to reach a certain accuracy is usually much lower [ ,55 J . TFOCS imple- 
ments such iterations and others like it but with various improvements. 

For large problems, e.g. images with a large number N of pixels, it is costly to hold the N x N opti- 
mization variable X in memory. To overcome this issue, our computational approach maintains a low-rank 
factorization of X. This is achieved by substituting the projection onto the semidefinite cone (the expensive 
step) with a proxy. Whereas V dumps the negative eigenvalues as in 

V{X) = ^max(Ai,0)niU*, 

i 

where Yli ^iUiU* (Ai > A2 > • • . > Aat) is any eigenvalue decomposition of X, our proxy only keeps the 
k largest eigenvalues in the expansion as in 



max(Ai, 0)iijii*. 



i<k 
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For small values of k — we use k between 10 and 20 — this can be efficiently computed since we only need 
to compute the top eigenvectors of a low-rank matrix at each step. Although this approximation gives good 
empirical results, convergence is no longer guaranteed. 



4.2 Setup 

To measure performance, we will use the mean-square error (MSB). However, since a solution xq is only 
unique up to global phase, it makes no sense to compute the squared distance between xq and the recovery 
xq. Rather, we compute the distance to the solution space, i.e. we are interested in the relative MSB defined 

as 

IICX0-X0II2 

mm r-;r . 

c:\c\ = l \\xo\r 

This is the definition we will adopt throughout the paper the Signal-to-Noise Ratio (SNR) is analogous, 
namely, SNR = logio(rel. MSB). 

Although our algorithm favors low-rank solutions, it is not guaranteed to find a rank-one solution. There- 
fore, if our optimal solution Xq does not have exactly rank one, we extract the rank-one approximation xqXq 
where xo is an eigenvector associated with the largest eigenvalue. We use a scaling such that ||xo|P = ll^^olP- 
Note that the £2 norm of the true solution is generally known since by Parseval's theorem, the £2 norm of 
Fxq is equal to ||xo||. Hence, observing the diffraction pattern of the object xq reveals its squared £2 norm. 



4.3 1-D simulations 

Phase retrieval for one-dimensional signals arises in fiber optics p6l[36l[37| , speech recognition [57], but 
also in the determination of concentration profiles in diffraction imaging p2[ . We evaluate PhaseLift for 
noiseless and noisy data using a variety of different 'illuminations' and test signals. 



4.3.1 Noisefree measurements 



In the first set of experiments we demonstrate the recovery of two very different signals from noiseless data. 
Both test signals are of length n = 128. The first signal, shown in Figure [3j a)) is a linear combination of 
a few sinusoids and represents a typical transfer function one might encounter in optics. The second signal 
is a complex signal, with independent Gaussian complex entries (each entry is of the form a + ih where a 
and b are independent A/'(0, 1) variables) so that the real and imaginary parts are independent white noise 
sequences; the real part of the signal is shown in Figure [3jb). 

Four random binary masks are used to perform the structured illumination so that we measure \Ax\'^ in 
which 

'Wi' 



A = F 



W2 
W3 



where each Wi is diagonal with either or 1 on the diagonal, resulting in a total of 512 intensity measure- 
ments. We work with the objective functional — ^(X)||2 + ATr(X) and the constraint X ^ Oto recover 
the signal, in which we use a small value for A such as 0.05 since we are dealing with noisefree data. We 



apply the reweighting scheme as discussed in Section 2.4 (To achieve perfect reconstruction, one would 



' Alternatively, we could use \\xoXq — a:oXo||F/|la;o3:^ol|f'' which gives very similar values. 
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(a) Smooth signal and its reconstruction (b) Random signal and its reconstruction (real part only) 

Figure 3: Two test signals and their reconstructions. The recovered signals are essentially indistin- 
guishable from the originals. 



have to let A — )• as the iteration count increases.) The algorithm is terminated when the relative residual 
error is less than a fixed tolerance, namely, ||^(xoXo) — &II2 < 10^^||6||2, where xq is the reconstructed 
signal just as before. The original and recovered signals are plotted in Figure [3f a) and (b). The SNR is 
87.3dB in the first case and 90.5dB in the second. 

We have repeated these experiments with the same test signals and the same algorithm, but using Gaus- 
sian masks instead of binary masks. In other words, the Wj's have Gaussian entries on the diagonal. It 
turns out that in this case, three illuminations — instead of four — were sufficient to obtain similar perfor- 
mance. This seems to be empirical support for a long-standing conjecture in quantum mechanics due to 
Ron Wright (see e.g. the concluding section of {bA}). The conjecture is that there exist three unitary op- 
erators Ui,U2, U3 such that the phase of the (finite-dimensional) signal x is uniquely determined by the 
measurements |?7i2;|, |f72x|, |?73x|. Our simulations suggest that one can choose Ui = F, U2 = FW, and 
C/3 = FW', where W, W are diagonal matrices with i.i.d. complex normal random variables as diagonal 
entries. The choice Ui = F, U2 = I, = FW was equally successful in our experiments, and is a bit 
closer to the quantum mechanical setting. Furthermore, we point out that no reweighting was needed, when 
we used six or more Gaussian masks. Expressed differently, plain trace-norm minimization succeeds with 
6n or more intensity measurements of this kind. 



4.3.2 Noisy measurements 

In the next set of experiments, we consider the case when the measurements are contaminated with Poisson 
noise. The test signal is again a complex random signal as above. Eight illuminations with binary masks are 
used. We add random Poisson noise to the measurements for five different SNR levels, ranging from about 
16dB to about 52dB. Since the solution is known, we have calculated reconstructions for various values of 
the parameter A balancing the negative log-likelihood and the trace norm, and report results for that A giving 



the lowest MSE. We implemented this strategy via the standard Golden Section Search 1 39 1. In practice one 
would have to find the best A via a strategy like cross validation (CV) or generalized cross validation (GCV). 
For each SNR level we repeated the experiment ten times with different random noise and different binary 
masks. 
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Figure 4: SNR versus relative MSE on a dB-scale for different numbers of illuminations with binary 
masks. The linear relationship between SNR and MSE is apparent. 



Figure |4] shows the average relative MSE in dB (the values of 10 logxo(rel. MSE) are plotted) versus the 
SNR. The error curves show clearly the desirable linear behavior between SNR and MSE with respect to the 
log-log scale. The performance degrades very gracefully with decreasing SNR. Furthermore, the difference 
of about 5dB between the error curve associated with four measurement and the error curve associated with 
eight measurements corresponds to an almost twofold error reduction, which is about as much improvement 
as one can hope to gain by doubling the number of measurements. 



4.3.3 Multiple measurements via oversampling 

It is well-known that for one-dimensional signals, oversampling does not result in unique solutions to the 



phase problem |44|. This might mainly be a theoretical issue without real practical consequences. Our 
numerical simulations confirm that this is not the case and demonstrate very clearly that oversampling is not 
useful for most one-dimensional signals. We apply one of Fienup's algorithms as discussed in Section 4.A 
in Q, and PhaseLift with reweighting to real-valued non-negative random signals of length 128. We use 
oversampling rates ranging from 2 to 5, and stop the algorithms when the relative residual error is less than 



Oversampling 


2 


3 


4 


5 


Relative MSE (Fienup) 


0.6811 


0.6548 


0.6180 


0.6150 


Relative MSE (PhaseLift) 


0.5451 


0.5343 


0.5299 


0.4930 



Table 1: MSE obtained by Fienup's algorithm and by PhaseLift from oversampled DFT measurements. 
Both methods produce guesses that fit the data and at the same time, are far from the true solution. This 
indicates that oversampling the DFT results in an ill-posed problem since we have several distinct 
solutions (probably infinitely many). 
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The results are displayed in Table [T] As can be seen, both algorithms find nearly perfect fits to the data 
and yet, they are very far from the original solution. Hence, we have a problem with multiple solutions. 
The average SNR of the "reconstructions" obtained via Fienup's method is around 4dB, which is just barely 
better than if we had used a random guess as a solution. The average SNR for the solutions computed by 
PhaseLift is only marginally better. 



4.4 2-D simulations 

We consider a stylized version of a setup one encounters in X-ray crystallography or diffraction imaging. 
The test image, shown in Figure [5ja) (magnitude), is a complex- valued image of size 256 x 256, whose 
pixel values correspond to the complex transmission coefficients of a collection of gold balls embedded in a 
medium. 



4.4.1 Noisefree measurements 

In the first experiment, we demonstrate the recovery of the image shown in Figure [5fa) from noiseless 
measurements. We consider two different types of illuminations. The first uses Gaussian random masks in 
which the coefficients on the diagonal of Wk are independent real- valued standard normal variables. We use 
three illuminations, one being constant, i.e. Wi = I, and the other two Gaussian. Again, we choose a small 
value of A set to 0.05 in 1 ^(X) |||+A Tr(X) since we have no noise, and stop the reweighting iterations 
as soon as the residual error obeys ||^(xoXo) — &II2 < 10^^||6||2- The reconstruction, shown in Figurejsjb), 
is visually indistinguishable from the original. Since the original image and the reconstruction are complex- 
valued, we only display the absolute value of each image throughout this and the next subsection. 

Gaussian random masks may not be realizable in practice. Our second example uses simple random 
binary masks, where the entries are either or 1 with equal probability. In this case, a larger number of 
illuminations as well as a larger number of reweighting steps are required to achieve a reconstruction of 
comparable quality. The result for eight binary illuminations is shown in Figure |5]^c). 



4.4.2 Noisy measurements 

In the second set of experiments we consider the same test image as before, but now with noisy measure- 
ments. In the first experiment the SNR is 20dB, in the second experiment the SNR is 60dB. We use 32 
Gaussian random masks in each case. The resulting reconstructions are depicted in Figure [6]; a) (20dB case) 
and Figure|6jb) (60dB case). The SNR in the 20dB case is 11.83dB. While the reconstructed image appears 
slightly more "fuzzy" than the original image, all features of the image are clearly visible. In the 60dB case 
the SNR is 47.96dB, and the reconstruction is virtually indistinguishable from the original image. 



4.4.3 Multiple measurements via oversampling 

Oversampling of two-dimensional signals is widely used to overcome the nonuniqueness of the phase re- 
trieval problem. We now explore whether this approach is viable. 

Here, we consider signals with real, non-negative values as test images, a case frequently considered in 
the literature, see e.g. [50-j52|. These images are of size 128 x 128. We take noiseless measurements and 
apply PhaseLift as well as Fienup's iterative algorithms [3, Section 4.A]. For each method, we terminate 
the iterations if the relative residual error is less than 10"'^ or if the relative error between two successive 
iterates is less than 10~^. 
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Figure 6: Reconstructions from noisy data via PhaseLift using 32 Gaussian random masks. 



• The simulations show that PhaseLift yields reconstructions that fit the measured data well, yielding a 
small relative residual error ||^(X) — y||2/||2/||2, yet the reconstructions are far away from the true 
signal. This behavior is indicative of an ill-conditioned problem. 

• The iterates of Fienup's algorithm stagnated most of the time without converging to a solution. At 
other times it did yield reconstructions that fit the measured data well, but in either case the recon- 
struction was always very different from the true signal. Moreover, the reconstructions vary widely 
depending on the initial (random) guess. 

Table [2] displays the results of PhaseLift as well as one of Fienup's algorithms, namely, the Error Reduction 
Algorithm [3, Section 4. A] (the other versions yield comparable results). The setup is this: we oversample 
the signal in each dimension by a factor of r, where r = 2,3, 4, 5. For each oversampling rate, we run ten 
experiments using a different test signal each time. The table shows the average residual errors over ten runs 
as well as the average relative MSE. The ill-posedness of the problem is evident from the disconnect between 
small residual error and large reconstruction error; that is to say, we fit the data very well and yet observe a 
large reconstruction error. Thus, in stark contrast to what is widely believed, our simulations indicate that 
oversampling by itself is not a viable strategy for phase retrieval even for non-negative, real-valued images. 



5 Discussion 

This paper introduces a novel framework for phase retrieval, combining multiple illuminations with tools 
from convex optimization, which been shown to work very well in practice and bears great potential. Having 
said this, our work also calls for improved theory, improved algorithms and a physical implementation of 
these ideas. Regarding this last point, it would be interesting to design physical experiments to test our 
methodology on real data, and we hope to report on this in a future publication. For now, we would like to 
bring up important open problems. 
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Algorithm | Oversampling 


2 


3 


4 


5 


P(X)-y||2/||y||2(Fienup) 


0.0650 


0.0607 


0.0541 


0.0713 


Relative MSB (Fienup) 


0.6931 


0.6882 


0.6736 


0.6878 


P(X)-y||2/||y||2(PhaseLift) 


0.0051 


0.0055 


0.0056 


0.0052 


Relative MSB (PhaseLift) 


0.4932 


0.4893 


0.4960 


0.4981 



Table 2: MSE obtained by Fienup's algorithm and by PhaseLift with reweighting from oversampled 
DFT measurements taken on 2D real-valued and positive test images. Fienup's algorithm does not 
always find a signal consistent with the data as well as the support constraint. (After the projection 
step in the spatial domain, the current guess does not always match the measurement in Fourier space. 
After 'projection' in Fourier space, the signal is not the Fourier transform of a signal obeying the 
spatial constraints.) Our approach always finds signals matching measured data very well, and yet the 
reconstructions achieve a large reconstruction error. This indicates severe ill-posedness since we have 
several distinct solutions providing an excellent fit to the measured data. 

At the theoretical level, we need to understand for which families of physically implementable structured 
illuminations does the trace-norm heuristic succeed? How many diffraction patterns are provably sufficient 
for our convex programming approach to work? Also, we have shown that our approach is robust to noise in 
the sense that the performance degrades very gracefully as the SNR decreases. Can this be made rigorous? 
Can we prove that our proposed framework is indeed robust to noise? Here, it is very likely that the tools 
and ideas developed in the theories of compressed sensing and of matrix completion will play a key role in 
addressing these fundamental issues. 

At the algorithmic level, we need to address the fact that the Ufting creates optimization problems of 
potentially enormous sizes. A tantalizing prospect is whether or not it is possible to use knowledge that the 
solution has low-rank, e. g. rank one, to design algorithms which do not need to assemble or store very large 
matrices. If so, how can this be done? Here, randomized algorithms holding up a sketch of the full matrix 
may prove very helpful. 

Finally, we would like to close by returning to another finding of this paper. Namely, oversampling the 

Fourier transform — this is the same as assuming finite support of the specimen — appears extremely prob- 
lematic in practice, even for real-valued nonnegative signals. To be sure, we have demonstrated that there 
typically exist very distinct 2D signals whose modulus of the Fourier transform nearly coincide, whatever 
the degree of oversampUng. In light of this extreme ill posedness, we have trouble understanding why this 
technique is used so heavily when it does not produce useful results in the absence of very specific a priori 
information about the image. Moreover, our concern is compounded by the additional fact that algorithms 
in common use tend to return solutions that depend on an initial guess so that different runs return widely 
different solutions. In the spirit of reproducible research, this calls for documented results making pubUcly 
available both data sets and software so that researchers can reproduce published results or results yet to be 
published. Of course, one might also be willing to rely on image priors far stronger than finite spatial extent, 
real-valuedness and positivity; they would, however, need to be specified. 
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